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Abstract 

We show that the decay of sinusoidal ripples on crystal surfaces, where mass transport is limited 
by the attachment and detachment of atoms at the step-edges, is remarkably different from the 
decay behavior that has been reported until now. Unlike the decreasing or at most constant rate 
of amplitude decay of sinusoidal profiles observed in earlier work, we find that the decay rate 
increases with decreasing amplitude in this kinetic regime. The rate of shape invariant amplitude 
relaxation is shown to be inversely proportional to both the square of the wavelength and the 
current amplitude. We have also carried out numerical simulations of the relaxation of realistic 
sputter ripples. 
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Nanoscale semiconductor devices hold potential for a wide range of applications in optics, 
electronics and information technology. The shape evolution of these structures at typical 
growth temperatures takes place through mass transport via surface diffusion. With increas- 
ing technological interest in smaller scale structures, a fundamental understanding of the 
kinetics of this phenomenon is important, since small feature sizes generally lead to more 
rapid evolution of the shape. Material parameters that govern atomic surface kinetics are 
efficiently determined from morphological relaxation of simple shapes, for example periodic 
sinusoidal surface ripples, which can be produced by lithographic pattering or are formed as 
a result of surface instabilities during growth, sputtering or etching. In the case of crystal 
surfaces above the roughening temperature and amorphous surfaces, the curvature based 
continuum theory developed by Herring and Mullins and the experiments supporting 
this theory for sinusoidal ripples, notably by Blakely and coworkers are considered classic 
studies in surface science. 

Surface evolution in the case of crystal surfaces below the roughening temperature takes 
place through diffusion of atoms on terraces and by their attachment and detachment from 
the step-edges. The surface energy in this case is determined by the formation and the 
interaction energies of steps and consequently acquires a cusp at facet orientations. Due to 
the presence of the strong non-linearity that arises from the cusp, extension of the classical 
theory to faceted surfaces is non-trivial. Earlier work on sinusoidal ripples, both continuum 
studies QflS and kinetic Monte Carlo (KMC) simulations fl Q focus on the regime where 
the kinetics of mass transport is determined by terrace-diffusion limited (TDL) kinetics. In 
the case of step attachment-detachment limited (ADL) kinetics, attention has been restricted 



to the relaxation of ID sinusoidal profiles 



where the line-tension of steps does not 



appear in the analysis. Many systems considered in recent experiments, notably 2D ripples 
on Si(OOl) Q, Cu(001)|ll|, Au(OOl) and Ni(OOl) Q, where the line-tension contribution 
should be significant due to the curvature of steps, are believed to be in the ADL regime. 

In this article, we consider the relaxation of 2D sinusoidal ripples (symmetric and elon- 
gated) and show that the decay rate of the amplitude increases with decreasing amplitude 
in the ADL case, for non-zero values of the formation energy (or line-tension) of steps. This 
is not an immediately obvious result, but can be established by numerical solution of the 
continuum surface evolution equations, scaling analysis and KMC simulations. We also 
carry out simulations of relaxation of realistic sputter ripples that consists of several Fourier 
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modes in the vicinity of a dominant peak in the surface spectrum. 

We start our continuum approach by relating the surface mass flux on the stepped surface 
to the gradient of the surface chemical potential /i, using the linear kinetic relation j = 
— C(Vh) V/U, where h denotes the height of the surface measured from a nominally flat 
reference plane, V denotes the surface gradient and C is the surface mobility parameter. 
The mobility parameter is taken to be of the form 

in which D s is the diffusion constant of adatoms, k is the effective rate of attachment 
and detachment of atoms at a step-edge, h s denotes the height of steps and C eq is the 
equilibrium adatom concentration ahead of a straight non-interacting step at temperature 
T. This functional form is obtained by taking the continuum limit of the terrace mass 
fluX from t he classic Burfc-CaWa-^ m o d e, for cWe surface st ep s fl Q. Tne 
relative importance of the terrace diffusion and the attachment-detachment processes can 
be obtained by comparing the effective length parameter I = D s /n with the terrace width, 
Ax = h s /\Wh\. For I <C Ax, the kinetics of mass transport is diffusion limited, whereas 
in the opposite limit, I ^> Ax, ADL kinetics is obtained. In the former case, the mobility 

D C 

parameter is simply a constant, D = £ Hp , while in the latter case, the mobility parameter is 
inversely proportional to the local slope. This recovers the result first derived by Nozieres 
for ADL kinetics. The expression in (JJJ describes a smooth transition between these two 
limiting cases. 

The driving force for surface diffusion can be obtained from the surface chemical potential, 
which can be written as 

where j(Vh) = 70 + A|V/i| + /5s|V/i| 3 is the usual form for the surface energy of a stepped 

surface in the small-slope limit. Here, 70 is the surface energy of the facet, /3i is related to the 

formation energy per unit length of a step and (3^ is related to the interaction between two 

steps. If the surface chemical potential is known, conservation of mass gives the equation 

for the evolution of surface height, 

dh 

- = -V ■ [C(Vh)V^] . (3) 

Because of the presence of a cusp in the surface energy, the chemical potential becomes sin- 
gular as \ Vh\ — > 0, which leads to difficulties in the numerical solution of © using standard 



discretization methods, unless the singularity is regularized. The solution, however can be 
accomplished without a regularization of the surface energy using a recent variational formu- 

n 

lation of surface evolution 14]. The basic idea of the formulation is to express the surface 
shape in terms of a Fourier series expansion and to use a variational principle to obtain 
coupled nonlinear ordinary differential equations (ODEs) for the expansion coefficients. If 
the initial surface shape is known, the evolution of these coefficients can be obtained using 
standard numerical integration techniques for ODEs. The convergence of the solutions can 
be assessed by varying the number of terms in the Fourier expansion. Below, we will analyze 
the relaxation of sinusoidal ripples in the case of ADL kinetics, first using the variational 
method, followed by a scaling analysis and KMC simulations. 

The influence of the step-edge barriers on surface relaxation can be conveniently discussed 
in terms of the dimensionless parameter e = Kh s /2D S , so that ADL kinetics would be 
observed as this parameter becomes small. The amplitude relaxation of a small amplitude 
2D sinusoidal ripple is plotted in Fig. 1 for three different values of this parameter, by 
choosing the ratio fli = fli/fls — 1. In all cases, we find that after an initial transient during 
which facets develop and grow at the extrema of the profile (the amplitude decreases by ~ 
20% in this period as shown in Fig. 1(a)), the shape relaxes in a self-similar manner at a rate 
that increases as the amplitude decreases in magnitude We have also studied amplitude 
relaxation for different values of the ratio f3\ and for sinusoidal profiles for which the ratios 
of the modulating wavelengths in the coordinate directions are different from unity. We 
find the decay behavior for finite values of (3\ to be fundamentally different from the decay 
for (3\ = 0. In the latter case, because of the absence of energetic cost for formation of 
steps, the decay shapes do not show any facets, but on the contrary, the top of the profiles 
become sharper in time, as shown in Fig. 1 (b). While the amplitude relaxation in this 
case is exponential (refer to Fig. 1(b)), any non-zero value of (3\ leads to a behavior similar 
to Fig. 1(a). These results indicate that the case ft 7^ is a singular perturbation of the 
special case with (3\ = 0. The observations regarding the functional form of relaxation do not 
change if the ratio of modulating wavelengths along the coordinate directions are changed 
- this result is similar to our earlier work in the TDL case fl4|. where decay behavior was 
shown not to depend on the aspect ratio of the ripples. We also note that ID sinusoidal 
profiles (like the 2D case in Fig. 1(b)) decay in an exponential manner when (3i = 0, which 
is in agreement with earlier work [a ll| where only the contribution from step interactions 
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was included in the surface chemical potential. 

Progress in understanding the surprising behavior in the case of ADL kinetics can be 
made using scaling arguments based on two key observations. First, for shape preserving 
relaxation, we can write the surface height as /i(x, t) = a(t)h(St), where x = k(xi, x%) denotes 



scaled coordinates and k = 2tt/X is the wavenumber of t 



re modulation considered in Fig. 1. 
our numerical simulations show 



Next, for (3\ ^ 0, just as in the case of TDL kinetics P, (| 
that the decay rate is insensitive to the value of the interaction strength fa during shape 
invariant relaxation. If the corresponding term in the surface chemical potential is assumed 
to be small a(t), and h(5£) satisfy 

v ; a(t) V %k B T J dxi \Vh\ d Xi x 3 dh Sj 

respectively, where the partial differential equation for h(x) does not depend on any physical 
parameters and the solutions are to be found in the domain < x~i,x 2 < 2tt. Since this 
equation is singular, it has to be solved by using its weak form or by employing variational 
methods introduced in Ref. Q|- the solution is very close to the faceted profile marked (iii) 
in Fig. 1(a). The ODE for a(t) can be readily integrated to obtain 



t , a(0) 2 / 2k B T \ 

a( t )=a(^l-- T , where T = -U- {j^-j . (5) 

The above analysis shows that the square of the amplitude should decrease linearly in time 
during shape invariant decay, which indeed is borne out by the numerical calculations using 
the variational approach (refer to the inset in Fig. 1(a)). A physical argument for the increase 
in decay rate with amplitude can be provided by noting that the material from the top of 
the profile that has to make its way to fill in the bottoms during relaxation, encounters 
fewer steps at later stages of decay compared to the early or intermediate stages. Since each 

n 

step presents a barrier to the atoms that are moving down [16], an increase in decay rate 
should then be expected with decreasing amplitude, or, alternatively, the total number of 
steps between the tops and bottoms of the profile. 

We now turn to an analysis of ripple relaxation using KMC simulations with a simple 
cubic solid-on-solid(SOS) model with near neighbor interactions. Here, we start with a 2D 
sinusoidal ripple shown in the inset in Fig. 2(a) and allow the atoms to perform hops at a 
rate determined by the number of neighbors and the barriers that they encounter during an 
atomic jump. The jump rate for an atom can be written as z/ exp[— (ne + E b )/k B T], where 



vq is an attempt frequency, n is the number of bonds attached to the atom in consideration, 
e is the bond strength and the barrier Eb is equal to the diffusion barrier Ed for hops on the 
same terrace, and is taken to be E^ + E s when an atom steps up or down a terrace, where 
E s is the additional barrier due to the Schwoebel effect |l6j |. Our calculations are similar 
to the work of Erlebacher and Aziz |7[ who considered the case E s = 0, but here we study 
surface relaxation as the parameter E s is varied. Since the attachment-detachment rate k 
in Eq. (0) decreases in an exponential manner with E s , the system would get closer to the 
ADL regime as this parameter increases in magnitude. 

The result of a typical relaxation simulation is shown in the inset in Fig. 2(a) for a system 
with an initial amplitude (ao) and wavelength of 6 and 80 lattice spacings, respectively. The 
profile shows a layer by layer decay, with very little step motion in other layers as the top 
and bottom layers are annihilated. The time dependence of the amplitude decay obtained by 
averaging twenty five statistically independent simulations at T = 0.41Tr, where T R = 0.62e 
is the roughening temperature, is shown in Fig. 2(a). When E s = 0, we recover the linear 
decay observed in earlier KMC simulations p, while the rate of decay increases with 
decreasing amplitude for larger values of E s . A fit in the case of E s /Ed = 1 shows that the 
amplitude decay is in agreement with the behavior predicted in Eq. Further evidence 
for the scaling given in Eq. (jSJ) can be found in Fig. 2(b), where for a fixed amplitude jl7|. 
the relaxation time scales with the square of the wavelength in the ADL, in contrast to cubic 
^obse^iotheTDLcaseSflQ. 

Finally, we consider the relaxation of realistic sputter ripples studied in recent experi- 
ments on Si and Cu surfaces jlfll . \v\ . While these ripples are characterized by a dominant 
wavelength, there is a spread in the surface roughness spectrum, due to the lack of perfect 
periodicity as shown in Fig. 3(i). Since the surface evolution equation Q is non-linear, 
a natural question to address is whether the decay of individual modes are influenced by 
the presence of neighboring modes. To this end, we have considered the evolution of the 
rippled surface in Fig. 3(ii) in the limit of ADL kinetics, by keeping track of approximately 
400 modes using the variational approach Q|. As the surface relaxes, smaller wavelength 
modes decay more rapidly, while the modes with larger wavelengths progressively dominate 
the spectrum. While a dominant wavelength can be identified at any given instant of time, 
its evolution neither follows the form given in (jSJ) nor shows an increasing rate of decay 
with amplitude. It can also be seen that some of the long wavelength modes, initially grow 
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in amplitude before they begin to decay. A similar trend has been recently observed on 



ripples on Cu surfaces 



11 1 . The above analysis shows that the non-linear coupling between 



modes on a typical rippled surface leads to complex relaxation behavior, and simple scaling 
rules for idealized sinusoidal corrugations cannot be directly applied to individual Fourier 
modes. We have also carried out similar analysis for TDL kinetics and found that the decay 
of dominant modes is not described by the linear relaxation behavior 0, 0, 7| of a single 
sinusoidal mode. It is possible that other metrics, for example the integrated amplitude of 
the power spectrum over all modes or the evolution of the peak and width of the spectrum 
might show trends that could be modeled with simple functional fits. We hope to explore 
these directions in the future. 

In conclusion, we have shown that the decay of sinusoidal profiles on crystal surfaces in the 
ADL limit, shows an increasing rate of relaxation with decreasing amplitude. This surprising 
behavior is attributed to kinetics of mass transport on stepped surfaces, in particular, the 
dependence of the surface mobility on the local density of steps. Analysis of rippled surfaces 
considered in experimental studies shows significant mode-coupling effects owing to the 
inherent non-linearity of the surface evolution equations. 
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FIG. 1: Amplitude relaxation for 2D profile h(x±,X2, t = 0) = ao cos(kx±) cos(kx2), where k is the 
wavenumber and ao = 0.01/fc. The relaxation is plotted as a function of t = DfykH, for (5\ = 1 
and in (a) and (b), respectively. As the parameter e = nh s /2D s in (a) decreases, the relaxation 
is expected to follow ADL kinetics. The inset graph shows the same curves, but with the vertical 
axis now set to a 2 (0, 0, t) /a§ in the region where the profile shows shape preserving relaxation - the 
curves are seen to be linear for small e, in agreement with Eq. Ij5]h . In contrast to (a), the decay 
in (b) is seen to be exponential in time, as evidenced by the linear dependence of the logarithm of 
the amplitude on time. The time sequences show facet formation (insets (ii) and (hi)) in (a), while 
the profile sharpens at the extrema in (b). 
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FIG. 2: (a) The normalized amplitude h(t)/ao plotted as a function of KMC time t for different 
ratios of the step barrier (E s ) to the diffusion barrier (Ed) for a 2D sinusoidal profile with amplitude 
and wavelength of 6 and 80 lattice spacings, respectively. In this calculation, we have chosen 
T = 0.41Tr and set Ed = e, the bond strength in the SOS model. The dashed lines for E s /Ed = 
and 1 represent fits to linear decay for TDL kinetics [j] and the decay law in Eq. ©l, respectively. 
The insets show the evolution of surface morphology for E s /Ed = 1, starting from the initial shape 
(i). (b) The normalized amplitude plotted as a function of KMC time for different wavelengths A 
(in lattice units) for the case E s /Ed = 1 at T = 0.41Tr. The inset shows that curves for different 
wavelengths nearly collapse to a single curve if time is scaled with the square of the wavelength. 
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FIG. 3: Scaled magnitudes of dominant Fourier modes |^4mn | = k\A mn \ in the expansion 
h(xi,X2,t) = Y^m n ^rn,n exp ik(mx\ + 71x2) for the shape of a rippled surface plotted as a function 
of rescaled time t = Df3^k A t for f}\ = 1 in the ADL regime. Inset (i) shows a typical AFM image of 
an unannealed sputter-rippled Cu surface. The size of the image is 5/jmx5ftm with an initial ripple 
wavelength and amplitude of ~ 300 nm and ~ 20 nm, respectively. Inset (ii) shows the unannealed 
sample of the same size used in the numerical study, with a dominant wavelength of 500 nm and 
an amplitude of 25 nm. 
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